/**
 * Name        : montecarlo.c
 * Author      : Edwin Luis Choquehuanca Mamani
 * Version     : 0.01
 * Copyright   : BSD
 * Description : Monte Carlo
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <assert.h>

#include <apr_general.h>
#include <apr_thread_proc.h>

#include <map>


#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 134456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349


/**
 * The number of running threads concurrently
 */
#define NUM_THREADS	64


typedef struct{
	int num;
	int ss;
}data;

static void* APR_THREAD_FUNC doit(apr_thread_t *thd, void *data);
double ran1(int *idum);
double gasdev(int *idum);

/**
 * Thread execution sample code
 * @remark Error checks omitted
 */
int main(int argc, const char *argv[])
{
	int N = 10000;
	int idum = 764;
	int a;
	double soma = 0.0;
	double temp;

	apr_status_t rv;
	apr_pool_t *mp;
	apr_thread_t *thd_arr[NUM_THREADS];
	apr_threadattr_t *thd_attr;
	int i;

	apr_initialize();
	apr_pool_create(&mp, NULL);

	/*for (a=0; a < N; a++) {
		temp = gasdev(&idum);
		soma += temp;
		printf("%f\n", temp);
	}*/

	/* The default thread attribute: detachable */
	apr_threadattr_create(&thd_attr, mp);
for (a=0; a < N; a++) {
	for (i = 0; i < NUM_THREADS; i++) {
		/* If the thread attribute is a default value, you can pass NULL to the second argument */
		rv = apr_thread_create(&thd_arr[i], thd_attr, doit, (void *)i, mp);
		assert(rv == APR_SUCCESS);
	}
}

	for (i = 0; i < NUM_THREADS; i++) {
		rv = apr_thread_join(&rv, thd_arr[i]);
		assert(rv == APR_SUCCESS);
	}

	apr_terminate();
	return 0;
}

/**
 * Thread entry point
 */
static void* APR_THREAD_FUNC doit(apr_thread_t *thd, void *data)
{
	int num = (int)data;
	gasdev(&num);
	//printf("doit: %d\n", num);
	//printf("rand: %d\n", gasdev(&num));
	apr_thread_exit(thd, APR_SUCCESS);

	return NULL;
}


double ran1(int *idum) {
	static long ix1, ix2, ix3;
	static float r[98];
	float temp;
	static int iff = 0;
	int j;
	if (*idum < 0 || iff == 0) {
		iff = 1;
		ix1 = (IC1-(*idum))%M1;
		ix1 = (IA1*ix1+IC1)%M1;
		ix2 = ix1%M2;
		ix1 = (IA1*ix1+IC1)%M1;
		ix3 = ix1%M3;
		for (j=1; j <=97; j++) {
			ix1 = (IA1*ix1+IC1)%M1;
			ix2 = (ix1*ix2+IC2)%M2;
			r[j] = (ix1+ix2*RM2)*RM1;
		}
		*idum = 1;
	}
	ix1 = (IA1*ix1+IC1)%M1;
	ix2 = (IA2*ix2+IC2)%M2;
	ix3 = (IA3*ix3+IC3)%M3;
	j = 1 + ((97*ix3)/M3);
	if (j > 97 || j < 1)
		;//erro

	temp=r[j];
	r[j] = (ix1+ix2*RM2)*RM1;
	return temp;
}

double gasdev(int *idum) {
	static int iset=0;
	static float gset;
	float fac, r, v1, v2;
	if (iset == 0) {
		do {
			v1 = 2.0*ran1(idum)-1.0;
			v2 = 2.0*ran1(idum)-1.0;
			r = v1*v1+v2*v2;
		} while (r >= 1.0);
		fac = sqrt(-2.0*log(r)/r);

		gset = v1*fac;
		iset = 1;
		return v2*fac;
	} else {
		iset = 0;
		return gset;
	}
}

